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ABSTRACT 

Estimates of the source sky location for gravitational wave signals are likely span areas ranging 
up to hundreds of square degrees or more, making it very challenging for most telescopes to search 
for counterpart signals in the electromagnetic spectrum. To boost the chance of successfully observ¬ 
ing such counterparts, we have developed an algorithm which optimizes the number of observing 
fields and their corresponding time allocations by maximizing the detection probability. As a proof- 
of-concept demonstration, we optimize follow-up observations targeting kilonovae using telescopes 
including CTIO-Dark Energy Camera, Subaru-HyperSuprimeCam, Pan-STARRS and Palomar Tran¬ 
sient Factory. We consider three simulated gravitational wave events with 90% credible error regions 
spanning areas from rsj 30 deg 2 to rsj 300deg 2 . Assuming a source at 200Mpc, we demonstrate that to 
obtain a maximum detection probability, there is an optimized number of fields for any particular 
event that a telescope should observe. To inform future telescope design studies, we present the maxi¬ 
mum detection probability and corresponding number of observing fields for a combination of limiting 
magnitudes and fields-of-view over a range of parameters. We show that for large gravitational wave 
error regions, telescope sensitivity rather than field-of-view, is the dominating factor in maximizing 
the detection probability. 


1. INTRODUCTION 

The detection of gravitational waves (GWs) from the 
inspiral and merger of binary black-hole systems (Ab¬ 
bott et al. 2016a, b) has marked the beginning of the 
GW astronomy era. With the advanced interferomet¬ 
ric detectors, GWs are expected to be observed from a 
number of additional source types including the mergers 
of binary neutron stars (BNS), and neutron star-black 
hole (NSBH) systems. For these other sources the pres¬ 
ence of matter from the neutron stars (NS) components, 
in such a violent interaction with its companion, makes 
it likely that electromagnetic (EM) emission will be gen¬ 
erated. The joint detection of a GW signal and its EM 
counterpart is therefore of particular interest. A coinci¬ 
dent EM observation together with the detection of GWs 
from these sources would lead to a deeper and more com¬ 
prehensive understanding of these objects. A successful 
EM follow-up observation triggered by a GW event can 
greatly improve the identification and localization of the 
host galaxy, bring us a more accurate estimate of the 
distance and energy involved, and provide better under¬ 
standing of the merger’s hydrodynamics and the local en¬ 
vironment of the progenitor (Bloom et al. 2009; Kulkarni 
& Kasliwal 2009; Phinney 2009). Additionally, knowl¬ 
edge obtained by EM follow-up observations may break 
the modeling degeneracies of binary properties, and con¬ 
firm the association between the GW source and its EM 
counterpart. Also, successfully locating the EM counter¬ 
part of a GW event could increase the confidence in the 
GW detection (Blackburn et al. 2015). 
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Among many other potential EM counterparts of GW 
events, it has been argued that kilonovae are a strong 
candidate for GW signals from BNS and NSBH merg¬ 
ers (Li & Paczynski 1998; Metzger & Berger 2011; Tan- 
vir et al. 2013; Li et al. 2016). We therefore focus 
our discussion on kilonovae only. Kilonovae are pre¬ 
dicted to produce an optical/infrared and isotropic quasi- 
thermal transient, and thought to originate from the hot 
neutron-rich matter ejected from BNS mergers. This 
ejection triggers r-process nucleosynthesis, the radioac¬ 
tive decay of which sustains the high temperature of 
the ejecta and powers the kilonovae. Kilonovae can 
be ^1000 times brighter than a nova (peak luminosity 
L ~ 10 40 erg/s) (Metzger et al. 2015), but given the ex¬ 
pected rate of BNS events, their corresponding distances 
make kilonovae relatively dim EM transients. In addi¬ 
tion, according to Metzger & Berger (2011), a kilonova 
maintains its peak luminosity only for hours to days af¬ 
ter merger. However, more recent calculations including 
the opacity of the r — process elements (Barnes & Kasen 
2013; Tanaka et al. 2014; Grossman et al. 2014) indi¬ 
cate that this timescale is likely to be days to weeks. 
Currently, there are three detections associated with the 
short-duration gamma-ray bursts GRB 130603B (Tanvir 
et al. 2013; Berger et al. 2013) and GRB 060614 (Yang 
et al. 2015; Jin et al. 2015) and GRB 050709 (Jin et al. 
2016). 

The sky location estimate for a GW event can cover a 
large fraction of the sky (>100 deg 2 (Singer et al. 2014)), 
posing a significant challenge to telescopes with fields of 
view (FOV’s) of order ~1 deg 2 trying to find counterpart 
signals in the EM spectrum. Observing and confidently 
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detecting kilonovae demands long exposure times even 
for powerful telescopes. It has been proposed that the in¬ 
frared sensitivity of the planned James Webb Space Tele¬ 
scope (JWST) may enable detections with short expo¬ 
sure times, thus compensating for its small FOV (Bartos 
et al. 2015b). The use of galaxy catalogs may also provide 
prior information on the direction of an event (Fan et al. 
2014; Bartos et al. 2015a; Gehrels et al. 2016). Nonethe¬ 
less, even with JWST and galaxy catalogs, the observa¬ 
tion time can still span an entire night, which may make 
target-of-opportunity BNS merger observations less at¬ 
tractive. 

Given that current and future EM telescopes will be 
subject to limited observational resources, we might 
ask how can we maximize the detection probability of 
a GW triggered kilonova signature. We have devel¬ 
oped an algorithm to answer this question, and here 
we present the results of a proof-of-concept demonstra¬ 
tion applied to four different telescopes including Subaru- 
HyperSuprimeCam (HSC) 1 (Miyazaki et al. 2012) , 
CTIO-Dark Energy Camera (DEC) (Bernstein et al. 
2012), Pan-Starrs 2 (Kaiser et al. 2002), and Palomar 
Transient Factory (PTF) (Law et al. 2009) for three dif¬ 
ferent simulated GW events. This algorithm takes as 
inputs GW sky localization information, and returns a 
guidance strategy for time allocation and telescope point¬ 
ing for a given EM telescope. In this paper we restrict 
the discussion to kilonovae from BNS mergers and leave 
the discussion of kilonova from NSBH for future study. 
This is in part due to the fact that our chosen GW data 
set from Singer et al. (2014) contains only BNS mergers. 
However, we remind the readers that there is no theoreti¬ 
cal restriction preventing us from applying our algorithm 
using other EM counterpart models and neither is it re¬ 
stricted to using only GW sky maps from compact binary 
coalescence (CBC) events. 

This paper is organized as follows: the methodology 
is introduced in Section 2 and we describe the imple¬ 
mentation in Section 3. The results obtained with the 
algorithm and a discussion of these results are given in 
Section 4. The possible future directions of this work are 
provided in Section 6, and our conclusions are presented 
in Section 7. 

2. METHODOLOGY 

As a proof-of-concept of our general EM follow-up 
method, here we have chosen to focus on kilonovae coun¬ 
terpart signatures. We have also adopted some simplifi¬ 
cations which will be relaxed in future studies. We as¬ 
sume that the available observation time is short com¬ 
pared to the luminosity variation timescale of a kilonova. 
This is validated by the fact that kilonova’s luminosity 
variation timescale is estimated to be ^days to a week, 
which is longer than a reasonable continuous EM obser¬ 
vation. This approximation allows us to assume that a 
kilonova has constant luminosity during the observation 
period. We also only consider the use of R-band lumi¬ 
nosity information, while the method presented can be 
extended to other regions of the EM spectrum. In real¬ 
ity, identifying a target EM counterpart requires tracking 
the object’s light curves, leading to several observations 

1 http://www. subarutelescope.org/Observing/Instruments/HSC/ 
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of the same point in the sky over several days. How¬ 
ever, in this work we deal with the problem of detec¬ 
tion, rather than identification. We consider only single 
observations and calculate optimized pointing directions 
and durations for a constrained total observation length 
and constant luminosity. Subsequent observations for 
the purpose of identifying variation within a field can be 
achieved by repeating our proposed observing strategy 
at some later time. 

In general, GW events could be localized to >100 deg 2 , 
which is much larger than a typical EM telescope’s FOV 
(<1 deg 2 ). We therefore do not consider the telescope’s 
rotation around its own axis, hence throughout this work, 
reference to single telescope pointing implies a rectangu¬ 
lar Charge-coupled device (CCD) image with edges par¬ 
allel to lines of longitude and circles of latitude. 

We also assume that the prior information from a GW 
trigger can be approximated as having independent sky 
location and distance probability distributions. Gener¬ 
ally, this may not be the case, however, our mathematical 
treatment can be greatly simplified under this assump¬ 
tion. We note that the sky maps available for our chosen 
dataset (discussed in Sec. 3) naturally lend themselves 
to this approximation since no distance information was 
computed as part of the rapid sky localization study re¬ 
ported in Singer et al. (2012) 3 . The final assumption is 
that the EM telescope can see every direction regardless 
of the location of the Sun, Moon, and the horizon. The 
issue of optimization of EM follow-up observations under 
the time-critical constraints such as those imposed by a 
source dipping below the horizon is explored in Rana 
et al. (2016). 

2.1. Bayesian Framework 

We use Dem to denote the successful detection of an 
EM counterpart. The probability of this occurring de¬ 
pends on the size of the selected telescope’s FOV co, the 
observed sky locations(cq J), and the exposure time r. 
The posterior probability of successful detection is then 
given by 

/*oo 

P(Dem\u,t, I) = / dN p(N\u,a,6,T,I). (1) 

Jn * 

Here, / is prior information that includes the selected 
telescope’s parameters, such as its photon collecting area 
A, filter and CCD efficiency. For a particular observation 
the number of photons N collected by the telescope is 

25 — m 

N = lO^s - x (tA) (2) 

where m is the apparent magnitude of the observed 
source. The threshold count V* is the criterion for de¬ 
tection determined by the signal-to-noise ratio (SNR) 
threshold, background noise and the selected telescope’s 
sensitivity. The value of N* is given by Eq. 2 with input 
values of m and r corresponding to the selected tele¬ 
scope’s detection threshold (see Table 1). The constant 
10 10 in Eq. 2 is the number of photons per second at 
m = 0. In practice, the value of V* should also ac¬ 
count for the change in background light accumulated 

3 The rapid sky localization algorithm BAYESTAR has since 
been updated to also provide rapid distance posteriors (Singer et al. 
2016 ). 
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TABLE 1 

Telescope parameters 


Telescope 

Aperture 

(m) 

FOV 

(deg 2 ) 

Exposure 

(<0 

Sensitivity 
(5-cr mag 
in R-band) 

TV */A 

(m- 2 ) 

DEC 

4.0 

3.0 

50 

23.7 

162.0 

HSC 

8.2 

1.13 a 

30 

24.5 

46.5 

Pan-Starrs 

1.8 

7.0 

60 

22.0 

930.5 

PTF 

1.2 

7.0 

60 

20.6 

3378.3 

LSST 

6.7 

9.6 

15 

24.5 

23.26 


a The full HSC FOV is 1.77 deg 2 but ~ 20% is used for calibration 


purposes. 


for different choices of observation time r, however, for 
simplicity we ignore this effect in this work. Since the 
number of photons expected from a target EM counter¬ 
part depends on its absolute magnitude M, distance R 
from the telescope, and how likely that the GW event 
is located within the field being observed, Eq. 1 can be 
expanded such that 

f'OG 

P(D em \u,t,I)= dN 
J TV* 

f dadS p(N\M,R,r,I)p(M\I)p{a,6,R\I). (3) 

J OJ 

The quantity P(N\M, R, r, I) is the probability of receiv¬ 
ing N photons from a source, given its absolute mag¬ 
nitude M, distance T, and observation time r, and is 
described by a Poisson distribution. Since we assume 
that the prior distribution on the distance to the target 
EM counterpart is statistically independent of the prior 
distribution on its sky location, Eq. 3 can be written as 


I dM I 


dR x 


P(Dem\w,t,I) = Pg w(^) x Pem{t) 


where, 


Pg wM 
Pem(t) 


- J p(a, S\I)da dS , 

= J dM J dR j dN 

p(N\M,R,r,I)p(R\I)p(M\I). 


(4) 

(5a) 

(5b) 


It should be noted that the GW sky localization infor¬ 
mation used for this work has been marginalized over dis¬ 
tance, meaning that the GW information represents a 2D 
error region projected onto the sky. However, since the 
distance estimate inferred from a GW observation will 
typically be poorly constrained and can be well approx¬ 
imated by a Gaussian distribution (Veitch et al. 2015; 
Singer et al. 2016), we assume a Gaussian prior with 
mean = 200Mpc and standard deviation = 60Mpc for 
the distance. In principle any form of positional informa¬ 
tion can be incorporated into our analysis and therefore 
our method can be adapted to include more realistic GW 
distance information. It is also possible that further con¬ 
straints from galaxy catalogs can also be incorporated 
into our method (Fan et al. 2014). 

We assume the least informative prior on peak lumi¬ 
nosity such that p(L\I) oc L~ 2 . It then follows that the 
prior on peak magnitude is given by 


where we assume M has a prior range of (—13,-8) as 
defined by the peak magnitudes of the models in Barnes 
& Kasen (2013). 

The probability of EM detection as defined in Eq. 1 
considers only one observing field. Given the size of a 
GW sky localization error region, and the typical size of 
an EM telescope’s FOV, the number of fields needed to 
be considered is >1. If an error region enclosing 90% of 
the GW probability covers S deg 2 and the EM telescope 
FOV is re deg 2 , the maximum number n of fields 4 re¬ 
quired to cover the error region at 90% can be estimated 
as n<S/uj. 

One might assume that observing as many fields as 
possible is optimal but we will show that telescope time 
is better spent by observing k fields where k lies in the 
range [1, n\. This occurs when it becomes more beneficial 
to observe a particular field for longer than observing a 
new field. 

For a given total observation time T, we are free to 
choose which fields we observe and the observation time 
allocated to each field. We represent these quantities by 
the vectors and {r^} respectively where k is the 

total number of chosen fields. Maximizing the detection 
probability of a kilonova amounts to finding the values 
of these vectors and the value of k which maximizes: 

P(D EM \k) = P(D em \{lo <*>}, {r®}, I) 

k<n 

= J2p(D e MldVfV). (7) 

i=1 

The choices of {r^} are subject to the constraint that 

k 

kT 0 + J2 T i k) = A (8) 

i=1 

where To represents the time required to slew between 
telescope pointings and/or perform CCD readout, and 
is equal to max(slew time, CCD readout time). We treat 
To as independent of the angular distance between point¬ 
ings. 

The expression we have for kilonova detection proba¬ 
bility (Eq. 7) as a function of the number of observed 
fields depends on our choice of field location and obser¬ 
vation time within each field. Given a number of fields k 
we begin choosing fields with a greedy algorithm, which 
will be described in section 3. Once the k fields have 
been chosen, they are represented by {uj^} and Eq. 7 
is maximized over the parameter vector {r^} to obtain 
the optimal kilonova detection probability. This is then 
repeated for each k in the range [1, n\ to find the optimal 
number of observed fields k. 

3. IMPLEMENTATION 

In this section we describe the processes for applying 
the GW sky localization information and generating the 
optimal observing strategy. The flow chart in Fig. 1 is a 
visual representation of this process. 

As shown in Fig. 1, we require information regarding 
the sky position of the GW source which we obtain from 
the BAYESTAR algorithm (Singer & Price 2016). This 
algorithm outputs GW sky localization information using 


p{M\I) oc 10"^ 


(6) 


4 In practice, n can be slightly larger due to overlapping fields. 
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TABLE 2 

Simulated GW event parameters. 


Event ID 1 

SNR 

90% Region 

(deg 2 ) 

Chirp Mass 

(M 0 ) 

28700 

16.8 

302 

1.33 

19296 

24.3 

103 

1.28 

18694 

24.0 

28.2 

1.31 


al The Event ID is that given to the events used in (Singer et al. 
2014) 


Fig. 1.— The process for generating an optimized observing strat¬ 
egy. Our algorithm takes two inputs: the GW sky localization 
information and a set of telescope parameters. After integrating 
over the number of received photons N, the source distance R , and 
the source absolute magnitude M, the algorithm returns the prob¬ 
ability of detection of an EM kilonova signal Pem as a function 
of the field observation time r. In parallel, for each choice of the 
total number of observed fields k : the algorithm selects the fields 
using a greedy algorithm. Based on the enclosed GW probability 
within each field the corresponding optimized observation times 
are computed using a Lagrange Multiplier approach. The total 
EM detection probability is then output for each choice of k. 


solution for the values {r^} that maximizes the detec¬ 
tion probability given by Eq. 7. This is subject to the 
constraints defined in Eq. 8 and the value of k that re¬ 
turns the highest detection probability is identified as the 
optimal solution. The analysis therefore guides us as to 
which subset of fields should be observed with the se¬ 
lected telescope and how much time should be allocated 
to each of those selected fields given the total observation 
time constraint. 


a HEALPix 5 coverage of the sky. Each HEALPix point 
corresponds to a value of the GW probability and repre¬ 
sents an equal area of the sky. BAYESTAR can rapidly 
(in ~10s) generate event location information and has 
been shown in Singer et al. (2014) to closely match results 
from more computationally intensive off-line Bayesian in¬ 
ference methods (Veitch et al. 2015). The simulated GW 
events used in this work are from BNS systems and taken 
directly from the dataset used in Singer et al. (2014). 

In this work, we consider follow-up observations using 
four telescopes (HSC, DEC, Pan-Starrs, and PTF) for 
three simulated representative GW events (see Table 2) 
which are studied assuming three total observation times 
T = 2,4,6hrs. For any telescope, GW event and total 
observation time the first stage of our procedure is to cal¬ 
culate the maximum number of fields n required to cover 
the sky area enclosed by the 90% probability contour of 
the GW sky region. 

In order to identify the possible observing fields for 
a given GW event, the greedy algorithm identifies the 
least number of HEALPix locations on the sky whose 
sum of GW probability is equal to the desired confi¬ 
dence level (in our case 90%). Then, assuming that 
each of those points represents the center of an observ¬ 
ing field we compute the sum of GW probability from 
each HEALPix point whose center lies within each of 
those fields. The field returning the maximum sum of 
the GW probability among those fields will be the first 
field. Subsequent fields are found by the same proce¬ 
dure, with the HEALPix points in the previous fields 
ignored. The summed probability within each field is an 
accurate approximation to the quantity Pgw(^) as de¬ 
fined in Eq. 5a. The n selected fields are labeled in the 
order with which they are chosen and hence their label 
indicates their rank in terms of enclosed GW probability. 
Therefore the first k < n fields represent our optimized 
choice of the values of {u^} in Eq. 7. 

As shown in Eqns. 4, 5b and 7, the detection proba¬ 
bility achieved by observing k selected fields can be ex¬ 
pressed as the sum of the product of the EM and GW 
probabilities in each field. For each value of k in the 
range [1, n] we apply a Lagrange Multiplier to find the 

5 http://healpix.sourceforge.net 


4. RESULTS 

In this section we present the results of our algorithm 
using our four example telescopes applied to the follow 
up of three simulated GW events. These events are taken 
from the dataset used in (Singer et al. 2014) and are des¬ 
ignated with the IDs 28700, 19296 and 18694. The error 
regions for these events each cover ~300 deg 2 , ^100 deg 2 , 
and ~30 deg 2 respectively and details of these events are 
presented in Table 2. These events were chosen to repre¬ 
sent the potential variation in sky localization ability of 
a global advanced detector network. We highlight that 
the actual injected distances for our chosen events are 
51Mpc, 27Mpc and 12Mpc respectively, while for our 
analysis we have assumed a distance of 200 Mpc for each 
event. However, the validity of this work is not under¬ 
mined. The injected events were originally simulated 
assuming a 2 detector aLIGO First Observational Run 
configuration. Our analysis assumes a 2(3) detector de¬ 
sign sensitivity aLIGO configuration. The SNR and sky 
localization for an event at 200 Mpc in the latter con¬ 
figuration is comparable (to within factors of ~few) to 
events at a few tens Mpc for the former configuration. 

Figure 2 shows the optimized tiling of observing fields 
obtained using the greedy algorithm approach. For each 
telescope we show sky maps of the GW probability over¬ 
laid with the 90% coverage tiling choices for the three 
representative GW events. The FOV’s of the telescopes 
range from 1.13 deg 2 to 7.1 deg 2 and as such the required 
number and location of tilings differ accordingly. The 
largest and smallest number of observation tilings are 
230 and 7 for the largest GW error region (ID 28700) 
using HSC and for the smallest error region (ID 18694) 
using either Pan-Starrs or PTF respectively. 

Each of the event maps are the result of an analysis 
assuming only 2 GW detectors. Without a third detec¬ 
tor the sky location of an event is restricted to a thin 
band of locations consistent with the single time delay 
measurement between detectors. This degeneracy is par¬ 
tially broken with the inclusion of antenna response in¬ 
formation resulting in extended arc structures. The third 
event that we consider (ID 18694) has sufficient SNR and 
suitable orientation with respect to the detector network 
that even with 2 detectors the sky region is well localized 
and is only partially extended. We consider this event to 
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be approximately representative of sky maps obtained 
from a 3 detector network. We also note that given the 
imperfect duty factors of both the initial and advanced 
detectors it is highly likely that future detections will be 
made whilst one or more detectors are offline. We there¬ 
fore use the first 2 example events (ID 28700 and 19296) 
as simultaneously representative of such a 2-detector sce¬ 
nario and of the potential situation in which a third de¬ 
tector is significantly less sensitive than the other two. 

Figures 3, 4 and 5 display the results of the simulated 
EM follow-up observations for the events labeled 28700, 
19296 and 18694 respectively. From the top panels to 
the bottom panels, the assumed total observation times 
are 6 hrs, 4hrs and 2 hrs. The plots on the left of the fig¬ 
ure show the optimized detection probability P(Dem|&) 
as a function of the total observed number of fields k. 
The plots on the right of the figure display the optimal 
time allocations corresponding to the value of k returning 
the highest detection probability (indicated by a circular 
marker in the detection probability plots on the left of 
the figure). The indices refer to the labels assigned to 
the fields when being chosen by the greedy algorithm. 

We highlight the asymptotic behavior of the detection 
probability curves in Figs. 3, 4 and 5. This is due to a 
particular feature of our algorithm and is explained as fol¬ 
lows. As the number of fields are increased we approach 
an optimal value k* where to observe an addition field 
with any finite observation time would actually reduce 
the detection probability. This occurs where the gains 
from an additional field are outweighed by the losses in¬ 
curred by reducing the lengths of the observations of the 
other fields. In this case, for a given value of k above 
the optimal value the optimal choice is to allocate r = 0 
to all fields with index greater than k*. If we know that 
we will allocate no time to these fields then we also have 
no need to slew to them or to readout from the CCD. 
Hence the optimal time allocations and also detection 
probability for values of k > k* remains constant at the 
maximum value. 

As a comparison, we introduce a second observing 
strategy for time allocation in which all the fields are 
observed with equal time. We call this method the equal 
time strategy. For each value of k , subject to the total 
observation time constraints and slew/readout time, the 
time allocated to each field is given by T/k — T 0 . The 
values of P(Dem|&) obtained using this strategy are plot¬ 
ted as dashed lines in the plots on the left of Figs. 3, 4 
and 5. The time allocations corresponding to the peaks 
of the dashed lines are plotted as fiat lines (constant equal 
values) in the plots on the right of the figures. The re¬ 
sultant maximal probabilities and corresponding optimal 
number of fields for both time allocation strategies are 
given in Table 3 for all simulated events and total ob¬ 
servation times. We also indicate the relative gains in 
detection probability obtained using our fully optimized 
(Lagrange multiplier) approach relative to the equal time 
strategy. 

In Fig. 6 we provide an insight into the detection po¬ 
tential of future telescopes and those not included in our 
analysis. We have computed the detection probability 
P(DEM\k) and corresponding optimal number of fields 
k as a function of arbitrary FOV and telescope sensitiv¬ 
ity. We define this sensitivity via the quantity N*/A, 


the number of photons per m 2 required for detection 
(see Eq. 2). For this general case we consider only a 
6 hr total observation of each of the 3 simulated events. 
For reference we include the 4 telescopes already con¬ 
sidered plus the proposed Large Synoptic Survey Tele¬ 
scope (LSST) (Abell et al. 2009) plotted with points 
indicating their locations in the FOV, telescope sensi¬ 
tivity plane. The relevant parameters for all telescopes 
included are given in Table 1. 

5. DISCUSSION 

The behavior of the detection probability as a function 
of the number of observed fields (as shown in Figs. 3, 4, 
and 5) shows that the 2 time allocation strategies pro¬ 
duce similar detection probabilities. In all cases the opti¬ 
mized approach gives marginally greater probability. As 
listed in Table 3 we see that for the particular cases ex¬ 
amined in this work the fully optimal approach leads to a 
typical gain of a few percent in detection probability over 
the equal time strategy. The biggest gains of ~ 5% are 
obtained for the DEC telescope. These relatively modest 
gains suggests that spreading observation time equally is 
close to optimal at the optimal number of observing fields 
k*. We also note that the number of fields at which the 
peak detection probability is achieved is similar for both 
strategies but always marginally lower for the equal time 
approach. 

Both strategies clearly indicate that with a given tele¬ 
scope, there exits a number of fields at which the prob¬ 
ability of a successful EM follow-up is maximized for a 
given GW event and a fixed amount of total observation 
time. In other words, exploring more or fewer fields than 
necessary will result in a decrease in detection proba¬ 
bility. Although one would expect the probability of a 
successful EM follow-up first increases with the number 
of fields observed, a drop occurs if so many fields are ob¬ 
served that only observations of short exposure time in 
some or all of the observed fields are allowed, given that 
the observation time is fixed. This trade-off between ex¬ 
ploring new fields and achieving increased depth within 
fields is broadly consistent with Nissanke et al. (2012). 
Overall, for the same event, the peak detection probabil¬ 
ities increase as the total observation time increases. For 
the same telescope, the detection probability increases as 
the size of error region decreases. This is expected since 
smaller error regions mean that a telescope requires fewer 
fields to cover that region at a given confidence. By ex¬ 
amining the curves in Figs. 3, 4 and 5, it can be seen 
that an increase in the total observation time shifts the 
position of the peak to larger numbers of fields, but does 
not change the general shape of the function. 

The time allocations shown in the plots on the right 
of Figs. 3, 4 and 5 are those computed for the optimal 
number of fields. In all cases our optimal strategy en¬ 
sures that proportionally more time will be assigned to 
those fields containing the greater fraction of GW prob¬ 
ability (the lower index fields). In general the range in 
observation times per field spans ~ 1 order of magnitude 
for a given optimized observation. A surprising feature 
of these distributions is that the optimal time alloca¬ 
tions can therefore differ by factors of a few per field 
with respect to the equal time distribution. However, 
both distributions result in very similar detection prob¬ 
abilities. Recently, Coughlin & Stubbs (2016) have pro- 
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(a) 90% coverage tilings used for the HSC telescope. 



(b) 90% coverage tilings used for the DEC telescope. 



xlO " 5 

8 

6 

4 


2 


1 


xlO " 4 

3 

2.5 
2 

1.5 


Right Ascension 


(c) 90% coverage tilings used for the Pan-Starrs telescope. 




xlO 


-4 



Right Ascension 


3 

2.5 

2 


l 

0.5 


(d) 90% coverage tilings used for the PTF telescope. 

Fig. 2.— The optimized locations of the observing fields covering 90% of the GW probability for each telescope and for each of the 3 
simulated GW events. Each row of plots corresponds to the 4 different telescopes (a) HSC, (b) DEC, (c) Pan-Starrs, and (d) PTF. Within 
each row we show the optimized observing field locations for the ~300deg 2 event (ID 28700) on the left, the ~100deg 2 event (ID 19296) 
in the center, and the ~30deg 2 event (ID 18694) on the right. In each plot the GW sky error is shown as a shaded region with the color 
bar indicating the value of posterior probability density. 













































7 




(a) EM detection probability and optimized observation times for T = 6 hour total observation times. 



(b) EM detection probability and optimized observation times for T = 4 hour total observation times. 




(c) EM detection probability and optimized observation times for T — 2 hour total observation times. 

Fig. 3.— The results of simulated EM follow-up observations for the ~300deg 2 GW event (ID 28700). We show the optimized EM 
detection probability as a function of the number of observing fields (left) and the allocated observing times for the optimal number of 
fields (right). The subfigures (a), (b) and (c) show results from 3 different total observation times for 6 hrs, 4 hrs and 2 hrs, respectively. 
For each total observation time the 4 solid curves in each plot correspond to the optimal time allocation strategy applied to each of the 4 
telescopes. The dashed lines show results for the equal time strategy. The solid markers and the circles indicate the number of observing 
fields at which the maximum detection probability is achieved using the optimal time allocation strategy and the equal time strategy 
respectively. 




























































(a) EM detection probability and optimized observation times for T = 6 hour total observation times. 




(b) EM detection probability and optimized observation times for T = 4 hour total observation times. 




(c) EM detection probability and optimized observation times for T — 2 hour total observation times. 

Fig. 4.— The results of simulated EM follow-up observations for the ~100deg 2 GW event (ID 19296). We show the optimized EM 
detection probability as a function of the number of observing fields (left) and the allocated observing times for the optimal number of 
fields (right). The subfigures (a), (b) and (c) show results from 3 different total observation times for 6 hrs, 4 hrs and 2 hrs, respectively. 
For each total observation time the 4 solid curves in each plot correspond to the optimal time allocation strategy applied to each of the 4 
telescopes. The dashed lines show results for the equal time strategy. The solid markers and the circles indicate the number of observing 
fields at which the maximum detection probability is achieved using the optimal time allocation strategy and the equal time strategy 
respectively. 
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(a) EM detection probability and optimized observation times for T = 6 hour total observation times. 




(b) EM detection probability and optimized observation times for T = 4 hour total observation times. 




(c) EM detection probability and optimized observation times for T — 2 hour total observation times. 

Fig. 5.— The results of simulated EM follow-up observations for the ~30deg 2 GW event (ID 18694). We show the optimized EM 
detection probability as a function of the number of observing fields (left) and the allocated observing times for the optimal number of 
fields (right). The subfigures (a), (b) and (c) show results from 3 different total observation times for 6 hrs, 4 hrs and 2 hrs, respectively. 
For each total observation time the 4 solid curves in each plot correspond to the optimal time allocation strategy applied to each of the 4 
telescopes. The dashed lines show results for the equal time strategy. The solid markers and the circles indicate the number of observing 
fields at which the maximum detection probability is achieved using the optimal time allocation strategy and the equal time strategy 
respectively. 
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(a) The EM detection probability and optimized number of fields for a general telescope observing event 28700 (~300deg 2 ). 



(b) The EM detection probability and optimized number of fields for a general telescope observing event 19296 (~100deg 2 ). 
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(c) The EM detection probability and optimized number of fields for a general telescope observing event 18694 (~30deg 2 ). 

Fig. 6.— Contours of EM follow-up performance of kilonovae as a function of the size of telescope FOV and sensitivity assuming a 6 hr 
total observation. We show results for the 3 simulated GW events (a) event ID 28700, ~300deg 2 , (b) event ID 19296 ~100deg 2 , and event 
18694 ~30deg 2 . For each event we plot contours of equal detection probability (left) and the corresponding optimal number of observing 
fields (right). Overlayed for reference on all plots are the locations of the telescopes considered in this work (including the proposed LSST). 
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TABLE 3 

The EM detection probability using both the optimal and equal time strategies 


Telescope 

Event ID 

Strategy 

EM detection probability (optimal number of fields) 

6 hrs 4 hrs 2 hrs 

6 h 

Relative Gain 
4 h 

2 h 

HSC 

28700 

19296 

18694 

LM a 

ET b 

LM 

ET 

LM 

ET 

66.4% 

65.0% 

78.1% 

77.0% 

93.7% 

92.9% 

(226) 

(198) 

(106) 

(103) 

(47) 

(47) 

58.0% 

57.0% 

75.1% 

73.7% 

92.7% 

91.5% 

(167) 

(155) 

(100) 

(93) 

(45) 

(44) 

41.6% 

41.3% 

65.3% 

63.7% 

89.1% 

87.4% 

(94) 

(90) 

(76) 

(71) 

(37) 

(36) 

1.4% 

1.1% 

0.8% 

1.0% 

1.4% 

1.2% 

0.3% 

1.6% 

1.7% 

DEC 

28700 

19296 

18694 

LM 

ET 

LM 

ET 

LM 

ET 

49A% 

48.0% 

71.5% 

67.2% 

85.3% 

83.0% 

(69) 

(60) 

(50) 

(39) 

(16) 

(16) 

4177% 

40.7% 

64.7% 

60.9% 

81.7% 

78.7% 

(56) 

(51) 

(41) 

(32) 

(16) 

(16) 

28.1% 

27.6% 

52.6% 

51.5% 

73.0% 

69.1% 

(35) 

(34) 

(19) 

(17) 

(16) 

(13) 

1.4% 

4.3% 

2.3% 

1.0% 

3.8% 

3.0% 

0.5% 

1.1% 

3.9% 

Pan-Starrs 

28700 

19296 

18694 

LM 

ET 

LM 

ET 

LM 

ET 

346% 

33.4% 

57.4% 

55.9% 

77.2% 

73.2% 

120)— 

(17) 

(12) 

(11) 

(7) 

(6) 

2875% 

28.0% 

50.1% 

49.1% 

71.5% 

67.8% 

-JT4]— 

(13) 

(10) 

(10) 

(7) 

(5) 

19.0% 

18.8% 

36.2% 

35.8% 

59.7% 

58.3% 

(9) 

(9) 

(8) 

(7) 

(5) 

(3) 

1.2% 

1.5% 

4.0% 

0.5% 

1.0% 

3.7% 

0.2% 

0.4% 

1.4% 

PTF 

28700 

19296 

18694 

LM 

ET 

LM 

ET 

LM 

ET 

1777% 

17.5% 

34.1% 

33.8% 

56.9% 

56.1% 

(9) 

(8) 

(7) 

(7) 

(4) 

(3) 

1370% 

12.9% 

25.8% 

25.6% 

50.1% 

49.2% 

(6) 

(6) 

(6) 

(5) 

(3) 

(3) 

7.1% 

7.1% 

14.6% 

14.6% 

34.9% 

33.5% 

(3) 

(3) 

(3) 

(3) 

(3) 

(3) 

0.2% 

0.3% 

0.8% 

0.1% 

0.2% 

0.9% 

0.0% 

0.0% 

1.4% 


a LM: Lagrange multiplier 
b ET: Equal time strategy 
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duced an analytic result for the distribution of observing 
times under a number of simplifying assumptions includ¬ 
ing a uniform prior on peak luminosity. They recommend 
that the time spent per field should be proportional to 
the prior GW probability in each field to the 2/3 power. 
We find that our numerical results are broadly consis¬ 
tent with this result but highlight again that the relative 
gains in detection probability are quite insensitive to the 
time allocation distribution. 

Figure 6 shows our results from another perspective, 
namely the optimal performance of a any existing or fu¬ 
ture telescope of arbitrary sensitivity and FOV. As men¬ 
tioned in Sec. 2 our treatment of the detection threshold 
criterion N* is simplified and hence these results should 
be treated as illustrative rather than definitive. However, 
as might be expected, a telescope with poor sensitivity 
and small FOV will be unlikely to detect an EM counter¬ 
part unless the GW event is particularly well localized. 

That LSST should explore more fields than PTF for 
these GW events may seem slightly confusing at first 
glance as LSST’s FOV is larger than PTF’s. This is 
because LSST is far more sensitive than PTF and can 
therefore explore as many fields as needed to cover the 
entire error region. In comparison PTF has to spend 
a considerable fraction of its total observation time on 
each of its observed fields, which results in PTF being 
only able to observe a limited number of fields. 

In general, the EM detection probabilities achieved are 
more sensitive to the telescope parameters for events 
with smaller error regions. For example, imagine that 
a telescope with a 2 deg 2 FOV and a R-band limiting 
magnitude of 20.3 in 30 sec (N*/A = 2273) needed to 
raise the detection probability P(T>em|&) by 10%. For 
event 18694, where the size of the 90% credible region is 
^30 deg 2 , it could either increase its FOV by a factor of 
« 2, or reduce its limiting magnitude to « 21.0 in 30 s. 
However, if the 90% credible region is ^300deg 2 as for 
event 28700, it would have to either further increase its 
FOV by a factor of ~ 2.5 or enhance its limiting mag¬ 
nitude to ~ 21.5 in 30 sec to have the same factor of 
improvement. 

In addition, the fact that at high sensitivity the con¬ 
tours in Fig. 6 appear to be almost flat indicates that 
as the size of the GW error region becomes larger, the 
FOV of a telescope has negligible impact on the detec¬ 
tion probability P(Dem |cj). For the design of future 
EM telescopes performing follow-up observations of GW 
events there will likely be a trade-off between sensitiv¬ 
ity and FOV. This result implies that for GW triggers 
with relatively large error regions, sensitivity, rather than 
FOV, is the dominant factor determining EM follow-up 
success. However, we remind the reader that this result 
is based on particular choices of source and correspond¬ 
ing prior on the source luminosity (see Eq. 6). Different 
choices may have different impacts on the outcome of our 
method. 

6. FUTURE WORK 

This work considers source sky error regions based 
solely on the information obtained from low latency GW 
sky localization (Singer & Price 2016). One may also 
include galaxy catalogs to further constrain source lo¬ 
cations within our existing Bayesian approach. In this 
work, the GW source distance was assumed to be statis¬ 


tically independent of its sky location. In the future we 
plan to use more realistic distance information (Singer, 
L. P. et al. 2016) therefore enhancing the effectiveness of 
our follow-up optimization strategy. 

Moreover, since telescopes are distributed at various 
latitudes and longitudes on the Earth, different tele¬ 
scopes are able to see different parts of the sky at dif¬ 
ferent times. As has been studied by (Rana et al. 2016) 
we would include the effect of observation prioritiza¬ 
tion when considering the diurnal cycle and for instances 
where GW sky error regions may pass below the horizon 
during a follow-up observation. Depending on the type 
of telescope considered we would also incorporate factors 
such as the obscuration of the source by the Sun and/or 
moon. 

In this work, we ask the question of detecting kilo- 
novae rather than identifying and characterizing them. 
The task of source identification is more demanding and 
will require the ability to differentiate our desired sources 
from contaminating backgrounds such as SNe and M- 
dwarf flares. One way to accomplish this is to perform 
multiple observations of the same fields. In this case a 
tentative detection is followed by an observation of the 
candidate for a period of time until the source’s light 
curve allows it to be classify as a kilonova or contam¬ 
inating noise (Cowperthwaite & Berger 2015; Nissanke 
et al. 2012). As mentioned in Sec. 2, simply repeating 
our proposed observations would enable the identifica¬ 
tion of variable objects for deeper follow-ups. However, 
a more involved procedure could incorporate light-curve 
information into our strategy therefore jointly optimizing 
the pointings used in both the detection and identifica¬ 
tion stages. This more complex strategy is left for future 
implementation. 

7. CONCLUSION 

In summary, we have demonstrated a proof-of-concept 
method for quantifying and maximizing the probability 
of a successful EM follow-up of a candidate GW event. 
We have applied this method to kilonovae counterparts 
but we emphasize that this method is versatile and ap¬ 
plicable to any EM counterpart model. We have shown 
that there exists an optimal number of fields on which 
time should be spent, and that observing more or fewer 
fields will result in a decrease in the detection probabil¬ 
ity. This analysis has been based on the assumptions of a 
static telescope with unconstrained pointings, a kilonova 
source at constant peak luminosity and the independence 
of statistical uncertainty between the distance and the 
GW trigger sky location. 

Our approach takes as inputs the GW sky localization 
information, and the selected telescope’s characteristics. 
The method selects the observed field locations with a 
greedy algorithm and then uses Lagrange Multipliers to 
compute the time allocation for those fields based on 
maximizing the detection probability. We have tested 
the algorithm by optimizing the EM follow-up observa¬ 
tions of the HSC, DEC, Pan-Starrs, and PTF telescopes 
for three simulated GW events. By comparing the results 
of our methods with the results of equally dividing the 
observation time amongst the observed fields, we have 
shown that both strategies return similar results, with 
our method producing marginally larger detection prob¬ 
abilities. 
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In addition, we have provided estimates for the EM 
follow-up performance of a general telescope of arbitrary 
sensitivity and FOV. These results indicate that in terms 
of telescope design, the likelihood of its success in the 
follow-up of kilonova signals is approximately indepen¬ 
dent of the FOV for reasonably sensitive telescopes. 

To extend this work, it may be helpful to consider in¬ 
cluding constraints from galaxy catalogs on source lo¬ 
cation. Also, the assumptions that the kilonova lumi¬ 
nosity is constant during the observation period, more 
realistic treatment of the source distance, and consider¬ 
ation of the dynamics of the telescope with respect to 
the source should also be investigated. Finally, the in¬ 
clusion of multiple observations of the same fields should 


be implemented to help distinguish kilonovae from con¬ 
taminating sources. 
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